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Abstract 


This article considers the problem of analyzing associations between power spectra of multi¬ 
ple time series and cross-sectional outcomes when data are observed from multiple subjects. 
The motivating application comes from sleep medicine, where researchers are able to non- 
invasively record physiological time series signals during sleep. The frequency patterns of 
these signals, which can be quantihed through the power spectrum, contain interpretable 
information about biological processes. An important problem in sleep research is drawing 
connections between power spectra of time series signals and clinical characteristics; these 
connections are key to understanding biological pathways through which sleep affects, and 
can be treated to improve, health. Such analyses are challenging as they must overcome the 
complicated structure of a power spectrum from multiple time series as a complex positive- 
dehnite matrix-valued function. This article proposes a new approach to such analyses based 
on a tensor-product spline model of Cholesky components of outcome-dependent power spec¬ 
tra. The approach flexibly models power spectra as nonparametric functions of frequency 
and outcome while preserving geometric constraints. Formulated in a fully Bayesian frame¬ 
work, a Whittle likelihood based Markov chain Monte Carlo (MCMC) algorithm is developed 
for automated model htting and for conducting inference on associations between outcomes 
and spectral measures. The method is used to analyze data from a study of sleep in older 
adults and uncovers new insights into how stress and arousal are connected to the amount 
of time one spends in bed. 

KEY WORDS: Bayesian Analysis; Coherence; Heart Rate Variability; MCMC; Multivariate 
Time Series; Sleep; Smoothing Spline; Spectral Analysis; Tensor-Product ANOVA; Whittle 
Likelihood. 

1 Introduction 

Innovations in data collection and storage have led to an increase in the number of biomedical 
studies that record multiple time series signals and outcome measures in multiple subjects. 
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For many time series, including common signals such as blood pressure, heart rate and 
electroencephalography (EEG), frequency patterns that are quantihed through the power 
spectrum contain important information about biological processes. Consequently, studies 
whose goals are to understand how underlying biological mechanisms are connected to be¬ 
havioral and clinical outcomes often require an analysis of associations between outcomes 
and power spectra of multiple time series. 


Our motivating application comes from a sleep study whose goal is to better understand 
the pathways that connect sleep to health and functioning. In the study, heart rate vari¬ 
ability (HRV) is recorded in subjects during a night of sleep. HRV is measured through the 
series of elapsed times between consecutive heart beats, and its power spectrum provides 


indirect measures of psychological stress and physiological arousal (Hall et al., 2007). Upon 
awakening, subjects reported subjectively assessed sleep outcomes, such as the amount of 


time slept during the night, which are associated with many aspects of well-being (Buysse 


2014). Understanding the association between the power spectrum of HRV during different 


sleep periods (i.e. beginning, middle and end of the night) and self-reported sleep outcomes is 
essential to understanding how stress connects sleep to health and, consequently, for guiding 
the use of treatments of poor sleep. 


In the biomedical literature, a two-stage approach is typically used to analyze such data. 
In the hrst stage, power collapsed within pre-selected frequency bands is estimated individu¬ 
ally for each time series (Malik et al. , [1996 ; Hall et al., 2004| ). A power spectrum is a function 
of frequency; power collapsed within a frequency band is an integral of the power spectrum 
over a range of frequencies, which converts the functional parameter into a scalar measure. 
In the second stage, classical statistical methods, such as ANOVA and linear regression, are 
used to evaluate associations between these band-collapsed spectral measures and outcomes. 
Such an approach has three major drawbacks. First, it is highly dependent on the frequency- 
band collapsed measures selected and there exists a hot debate as to which measures should 


be considered and/or how they should be interpreted (Burr, 2007). Ideally, an analysis of 
such data should provide global measures that can be used to understand the entire system 
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while also providing a means to conduct inference on any frequency band-collapsed measure 
of potential interest. Second, estimated power is treated as if it were not an estimate but the 
true unknown parameter, leading to inaccurate inference. Finally, band-specific frequency 
measures are estimated for each time series separately, inhibiting the evaluation of patterns 
across series. For instance, in our motivating example, this two-stage approach does not 
examine how the coherence in HRV between the beginning and end of the night is connected 
to sleep outcomes. 


In the statistics literature, a considerable amount of research has been conducted on 


methods for analyzing functional variables, a thorough review of which is given by Wang 


et ah (2016). Included in this body of work are methods for analyzing associations between 


power spectra and outcomes when one time series is observed per subject (Stoffer et al. 


2010 Krafty and Hall, 2013). When one observes multiple time series per subject and 


interest lies in frequency patterns both within each series and across different series, the 
problem becomes considerably more challenging. This is the case in our motivating study, 
where we are interested not only in stress and arousal during particular periods of sleep, but 
also in their persistence and coherence across periods. While the power spectrum from a 
single time series is a positive real-valued function of frequency, the power spectrum from 
multiple time series is a positive-definite Hermitian matrix valued function of frequency. An 
analysis of associations between power spectra from multiple time series and study outcomes 
must be able to flexibly model associations while preserving this positive-definite Hermitian 
structure. 


Efficient nonparametric methods that preserve the positive-definite Hermitian structure 
of spectral matrices have been developed for the simpler, classical problem of estimating the 
power spectrum of a multivariate time series from a single subject by modeling Cholesky 


components of spectral matrices as functions of frequency (Dai and Guo, 2004; Rosen and 


Stoffer, 2007 Krafty and Collinge, 2013). In this article, we extend this framework to develop 


a new approach to analyzing data from multiple subjects that models Cholesky components 
as functions of both frequency and outcome. Rather than being curves as functions of 
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frequency, components of spectral matrices under the proposed model are surfaces. Changes 
in these surfaces with respect to the outcome provide nonparametric measures of association 
between outcomes and power spectra. The proposed method is formulated in a fully Bayesian 
framework; a MCMC algorithm based on the Whittle likelihood, or the asymptotic likelihood 
derived from the Fourier transform of the data, is developed for model htting and inference. 
The method allows one to evaluate the entire outcome-dependent power spectrum and to 
conduct nonparametric inference on the association between the outcome and any function 
of the power spectrum. 

The rest of the article is organized as follows. Our motivating application, the AgeWise 
Sleep Study, is discussed in Section A review of spectral analysis in the classical setting, 
where data are observed from a single subject, is given in Section]^ The proposed methodol¬ 
ogy for analyzing time series from multiple subjects is presented in Section]^ The proposed 
method is used to analyze data from the motivating application in Section and some hnal 
remarks are offered in Section El 


2 The AgeWise Sleep Study 


An estimated 43% of older adults report problems initiating or maintaining sleep (Foley 


et ah, 1995). Poor sleep in older adults has been linked to depression, heart disease, obesity. 


arthritis, diabetes and stroke (Foley et ah, 2004). With medical and scientihc advances 


leading to an increase in the world’s elderly population, the consequences of poor sleep 
in older adults pose a major public health concern. The AgeWise study is a NIH-funded 
Program Project conducted at the University of Pittsburgh that seeks a better understanding 
of causes, effects, and treatments of poor sleep in older adults. Towards this goal, we consider 
A^ = 108 men and women between 69-89 years of age who were observed during a night of 
in-home sleep. Two types of data were collected in each subject. First, subjects were 
observed during the night through ambulatory polysomnography (PSG), or the continuous 
collection of electrophysiological changes that occur during sleep. Second, upon awakening. 


subjects completed the Pittsburgh Sleep Diary (Monk et ah, 1994) to record self-reported 
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sleep outcomes during the night. 


As previously discussed, HRV is the series of elapsed times between heart beats. It is 
of interest to researchers, as it reflects neurological control of the heart, and through this 
capacity, its power spectrum provides indirect measures of stress reactivity and arousal. The 
PSG used in the study included an electrocardiograph (EGG) to monitor heart activity. The 
EGG was used to locate the timing of heart beats, which were then differenced, detrended, 
cubic spline interpolated, and resampled at 1 Hz to compute HRV series throughout the 
night. 


During the night, the body cycles through two types of sleep: rapid eye movement 
(REM) and non-rapid eye movement (NREM) sleep. In NREM sleep, which contains deep- 
sleep, the parasympathetic branch of the autonomic nervous system that is responsible for 
unconscious actions and stimulates the body to “rest-and-digest” dominates the sympathetic 
branch, which drives the “flight-or-£ght” response. Parasympathetic nervous system activity 
during NREM is hypothesized to be responsible for many of the rejuvenating properties of 


sleep (Siegel, 2005). However, physiological activity during NREM sleep is not constant; in 
good sleepers, the amount of parasympathetic activity during NREM increases throughout 


the night (Hall et ah, 2004). To enable an analysis that can evaluate autonomic nervous 
system activity and its changes during the night, we consider 3 HRV time series per subject 
(at the beginning, middle, and end of the night) by extracting the hrst 5 minutes of HRV 
from the hrst three periods of NREM sleep. Data from two subjects are displayed in Figure 

□ 


The goal of our analysis is to understand how the power spectrum of HRV over the 
three periods of NREM are connected to self-reported sleep. We focus on one particular 
self-reported sleep measure derived from the Pittsburgh Sleep Diary: time in bed (TIB). 
TIB is dehned as the elapsed time between attempted sleep and hnal wakening. It serves 
as an upper bound for the amount of time spent asleep during the night, which has been 
linked to heart disease, hypertension, impaired neurobehavioral performance and mortality 


(Buysse, 2014). The reported TIB from our sample has a mean of 477.99 minutes and a 
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Figure 1: Detrended HRV time series during the first three periods of NREM from two 
subjects. Subject 1 reported a TIB of 357.67 minutes and subject 2 reported a TIB of 
521.00 minutes. 

standard deviation of 71.32 minutes. The resulting data for analysis consist of three epochs 
of HRV, one during each of the hrst three periods of NREM sleep, and self-reported TIB 
from each subject. 


3 Methodological Background: Spectral Domain Anal¬ 
ysis 

Before introducing our proposed model for the spectral analysis of multiple time series from 
multiple subjects in Section in this section we present background on spectral analysis in 
the classical setting, where data are observed from a single subject, for both univariate and 
multivariate time series. 
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3.1 Univariate Time Series 


3.1.1 Population Parameters 

Spectral domain analysis focuses on the cyclic behavior of time series data. An alternate 
approach is time domain analysis wherein the relationship between the data at different time 
lags is the focus. For stationary time series, the main time domain tool is the covariance 
between a current value of the series, say X^, and the value of the series h time units prior, 
say Xt-h- The autocovariance is a function of lag, and is given by 

7(/i) = Cov(Xt, Xt-h), h = 0, ±1, ± 2 ,... . 

If 7 (h) is absolutely summable (which it is for ARMA models, for example), then there is a 
duality between the power spectrum, given by 

00 

/M= 7 (/i) exp {—27iiuh ), a; G M, 

h=—oo 

and the autocovariance function, namely, 

/. 1/2 

l{h) = / f{oj) exp {27riujh) doo, h = 0, ±1, ±2,..., (1) 

J-l/2 

as the inverse transform of the power spectrum. The relationship is the same as that of a 
characteristic function and a probability density. Consequently, the information contained 
in the power spectrum is equivalent to the information contained in the autocovariance 
function. If we are concerned with lagged behavior, working with 7 (/i) is more informative; 
if we are concerned with cyclic behavior, as is the case of HRV where cyclical behavior 
provides interpretable physiological information, working with f{u) is more informative. 

The power spectrum is nonnegative and we assume that it is positive, so that f{u) > 0 
for all u. In addition to being positive, / has two other restrictions as a function of frequency. 
By the nature of the Fourier transform, it is periodic such that f{u) = f{u + 2n), and it 
is a Hermitian function, or an even function, where f{u) = f{—u). Consequently, f{u) is 
usually displayed only for u G [0,1/2]. 



Putting /i = 0 in Q yields 

/•1/2 

7 ( 0 ) = Var(Xt) = / f{uj)duj, 

7-1/2 

which expresses the total variance of the time series as the integrated power spectrum. In 
particular, we may think of f{u) doo as the approximate variance in the data attributed to 
frequencies in a small band of width du around uj. It is common to view spectral analysis 
as an analysis of variance (ANOVA) of time series data with respect to frequency. In fact, 
the power spectrum is a density of variance rather than of probability. 

3.1.2 Estimation 


The nonparametric estimation of / from an epoch of length n, Xi,... ,X„, can begin by 
considering the discrete Fourier transform (DFT) 


Ym = n ^ Xt exp{-27iiumt), 


t=i 


where Um = m/n are the Fourier frequencies. When n is large, FA are approximately 
independent mean-zero complex normal random variables with variances f{ujm) for m = 
1,... ,M, M = [(n — 1)/2J, and Ym = Yn-m (Shumway and Stoffer 


2011 


Appendix C). 


Consequently, the periodogram \Ym\ provides approximately unbiased but noisy estimates 
of f{oJm)- Consistent estimates can be obtained by smoothing the periodogram across fre¬ 


quency using tools such as local averaging (Shumway and Stoffer, 2011, Chapter 4.5), splines 


Pawitan and O’Sullivan, 

1994 

), and wavelets ( 

Moulin 

1994 


Our estimation approach for multiple time series from multiple subjects, which we develop 
in Section is based on Bayesian splines. To motivate the development of the new method¬ 
ology, in this subsection we discuss first a Bayesian smoothing spline model for univariate 
spectral analysis, then discuss a low-rank approximation. Smoothing spline estimation bal¬ 
ances the fit of a function to observed data with a roughness-based measure of regularity. 
The Bayesian formulation of smoothing splines was first discussed in the case of Gaussian 


observations by Kimeldorf and Wahba (1970) and Wahba (1978). Under the Bayesian for¬ 


mulation, the likelihood provides a measure of fit to observed data, and regularity is imposed 
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through a mean-zero Gaussian prior on the functional parameter, which induces a prior for 
the roughness of the function. 


For spectrum estimation, the large sample distribution of fdn provides the Whittle like¬ 


lihood (Whittle, 1953, 1954) 


M 


m=l 

We adopt generic notation throughout this article where Y will denote all DFT data. Since 
/ is a positive function, log / is modeled rather than / itself to avoid constraints. Although 
general measures of regularity can be considered, we focus on measuring the roughness of a 
function through its integrated squared second derivative, 

^ 1/2 

'P (log /) = / {[log /]" (w) duj. 

Jo 

The specihcation of the prior distribution begins by decomposing log/ into a linear part 
(which is in the null space of V) and a nonlinear part. To dehne the prior distribution for 
the nonlinear part, consider the reproducing kernel of the seminorm dehned by V 

/•1/2 

J{ui, Uj) = {Ui- {Uj - z/)_^ du, 


where (u), = max (z/, 0) (Gu 


2013 


Section 2.3.1). For the Bayesian smoothing spline model. 


the prior distribution for the log-spectrum can be formulated as 

M 

. 


log / (w) = ai -F a 2 UJ + ^ J{u,u;j)zj 

i=i 


where z = {zi,...,zm)' ~ A^(0, is independent of a = (ai,ai)^ N{0,alh), J = 

{J{u}i, ojj)} is the M X M matrix of J evaluated at the Fourier frequencies and I 2 is the 2x2 
identity matrix. The reproducing property of the kernel J provides a simple form for the 


roughness of the log-spectrum, V (log/) = z'Jz (Gu, 2013, Ghapter 2), from which it can 
be seen that the prior distribution on the coefficients z induces a prior on the roughness of 
the spectrum where V (log/) ~ t^Xm^ cind Xm denotes a chi-squared random variable on 
M degrees of freedom. The smoothing parameter > 0 balances the smoothness of the 
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estimator to its fit to the data such that, conditional on the Bayes estimate (posterior 
median) approaches a linear function as —>■ 0 and interpolates as —)■ oo. 

Two approaches may be taken for estimation and inference under the Bayesian model: 
empirical Bayes and fully Bayes. In the empirical Bayes approach, a data driven method, such 
as generalized cross-validation (GCV) or generalized maximum likelihood (GML), is used to 
select the smoothing parameter r^. The log-spectrum is then estimated from the posterior 
distribution conditional on r^, and its median is equivalent to the frequentist smoothing 


spline obtained by minimizing a penalized Whittle likelihood as ct„ —)■ cx) (Gu, 1992 Qin 
and Wang, 2008). In the fully Bayesian approach, is treated as a random variable with 


a noninformative prior, and inference is conducted averaged over the posterior distribution 
of (Speckman and Sun, 2003 Grainiceanu et ah, 2005). Our proposed methodology will 


adopt the later approach, as discussed in Section 

The presented smoothing spline model contains a large number of coefficients, which 
can impede computation and limit practicality. Low rank approximations, such as those 


obtained by using a subset of the kernel functions J{-,ujj) (Gu and Kim, 2002) or another 


set of basis functions contained in the column space of J (Wood, 2006), can be used to 


ease computational burden without sacrihcing model £t. Here, we consider the basis formed 
from the scaled eigenvectors of J, which has been used for power spectrum estimation by 


Wood et ah (2002), Rosen et ah (2009) and Rosen et ah (2012). This basis can model 


smooth functions with a relatively few number of basis functions, provides a diagonal prior 
covariance structure that aids computation, and maintains the intuitive interpretation of the 


prior distribution regularizing roughness as measured through V (Nychka and Gummins 


1996). To formulate this low-rank approximation, we first consider an equivalent formulation 


of the Bayesian smoothing spline model for log / at the Fourier frequencies. Let log / = 
[log/(a;i),..., log/(a;M)]^ be the log-spectrum evaluated at the Fourier frequencies. Further, 
let J = VjDjVj be the spectral decomposition of J, Qj = VjD^p and c = D^j’^V'jZ. Then 
an equivalent model at the Fourier frequencies is 


log / = Lja + Qjc, 
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where Lj = | ; 1m is the M-vector of ones, u) = (cui,..., um)', and c rsj 


The eigenvectors contained in the columns of Vj are in increasing order of roughness 


and the eigenvalues contained in the diagonal matrix Dj decay rapidly (Eubank, 1999). 


Smooth functions can be accurately modeled through the first several nj < M to provide 
a low-rank approximation. With a slight abuse of notation to avoid the need to introduce 
further variables, we let Qj represent the M x nj matrix of the hrst rij eigenvectors with 
corresponding coefficients c ~ (0, 

The selection of nj provides a compromise between low-rank computational feasibility 
and loss of flexibility relative to full M-rank model. An intuitive measure of the loss of 
flexibility is the fraction of the total variance of the covariance matrix J explained by its 
first nj eigenvectors (EVE), or the sum of its first nj eigenvalues divided by the total sum 


of all of its M eigenvalues. Wood et ah (2002), Rosen et ah (2009) and Rosen et ah (2012) 


suggest using nj = 10 basis functions, which, in each of their settings considered, equates 
to an EVE of 97.975%. Our empirical findings support this suggestion, and we recommend 
selecting rij to achieve a 97.975% EVE. Under this rule, nj = 7 for n = [15,18], nj = S for 
n = [19,22], nj = 9 for n G [23,40] and rij = 10 for n G [41,10^]. 


3.2 Multivariate Time Series 

3.2.1 Population Parameters 

The ideas presented in the univariate case generalize to the multivariate case wherein we 
observe a P-dimensional vector-valued time series, say {Xj}. Under stationarity, the auto¬ 
covariance function is a P x P matrix given by 

T{h) = Cov{Xt,Xt-h), h = 0, ±1, ±2,... . 

If l|r(h)|| < oo, the spectral density matrix of the series Xt is given by 

OO 

/M= r(h) exp (—27ria;h), a; G M 

h=—oo 
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and the inverse relationship is 

/• 1/2 

r(^) = / f{^) exp{27iiujh) du, h = 0, ±1, ±2,... . 

J-l/2 

For each cj G M, /(w) is a P x P non-negative definite Hermitian matrix with the diagonal 
elements, fpp{oj) ior p = 1,..., P, being the spectra of the individnal components, and the 
off-diagonal elements, fqpioj) ioi p = 1,..., P, being the cross-spectra. Thronghont this 
article, we assnme that /(w) is non-singnlar for all a; G M. As in the nnivariate case, / is a 
periodic and Hermitian fnnction of freqnency were, for matrix-valned fnnctions, Hermitian 
as a fnnction of freqnency is defined as fioj) = f*{—u), and /*(a;) is the complex conjngate 
of f{uj). 


An important example of the application of the cross-spectrnm is to the problem of 
linearly predicting one of the component series, say Xqt, from another component, say Xpt. 
A measnre of the strength of snch a relationship is the sqnared coherence fnnction defined 
as 


Plpi^) = 




fqq{^)fpp{^) 

This is analogons to conventional sqnared correlation between two finite-variance random 
variables; e.g., 0 < Pqp{uj) < 1. This analogy motivates the interpretation of sqnared coher¬ 
ence as the sqnared correlation between two time series at freqnency u. These ideas extend 
in an obvious way to the concept of multiple coherence and partial coherence functions ob¬ 
tained from the full spectral matrix in much the same way that multiple correlation and 
partial correlation can be obtained from a covariance matrix. Full details of these results 


may be found in Shumway and Stoffer (2011, Chapters 4 & 7). 


3.2.2 Estimation 


In the multivariate setting, let 

n 

Ym = Xt exp{-27iiUrnt) 

t=l 
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be the DFTs of the data. In this case, the Whittle likelihood is 

M 

L{Y I /) ^ l[\r\uJ\exp{-YU-\ujjYm}, 

m=l 

and the periodogram YyrXm is an approximately unbiased but noisy estimate of f{oJm), from 
which consistent estimates can be obtained by smoothing. 


While in the univariate setting the spectrum is smoothed on the logarithmic scale to 
preserve positivity, Cholesky components of spectral matrices can be smoothed to preserve 


positive-dehniteness in the multivariate setting (Dai and Guo, 2004 Rosen and Stoffer, 2007 


Krafty and Collinge, 2013). The modihed Cholesky decomposition assures that, for a spectral 


matrix f{u), there exists a unique P x P lower triangular complex matrix 0(a;) with ones 
on the diagonal and a unique P x P positive diagonal matrix 4/(a;) such that 


= 0(a;)T-Xa;)0*(a;). 


There are P^-Cholesky components to estimate: 3? {Qki} and {0^} ior k > i = 1,..., P — 
1, and for A: = 1,..., P. Since the diagonal terms > 0, we model log Letting 

Oki = [0fcr(a;i),..., Okii^Ju)]' and log^^^ = [log ..., log we model: 


^{0ki} = 

Ljdrki Qj^ki") 

k>i = l. 

,...,P-1 

(2) 

^ {0ke} = 

LjOfiki Q J^iki") 

k>i=l, 

...,P-1 

(3) 

logV'fcfc = 

LjCldkk Q J^dkk") 

k = 1, . . . 

,P, 

(4) 

where Crki ~ Cik£ 

~ Cdkk 

~ N{0,T^kkQ, (^rke ~ 

A^(0, 0 -^/ 2 ), aike ~ 


A^(0, alh) and aakk iV(0, (j^/ 2 ). Throughout this article, r, i and d are used to denote 
coefficients for real components of 0, imaginary components of 0 and the logarithm of the 
diagonal components of 4/“^, respectively. 


4 Methodology: Replicated Multiple Time Series 

The primary question considered in this article is how to assess the association between the 
power spectrum of P-variate time series of length n, {Xji,... ,Xjn}, and real-valued static 
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variables, Uj, observed from j = 1,... ,N independent subjects. In the motivating study, 
there are N = 108 participants, Uj is self-reported TIB, and Xjt are time series of HRV 
during the first P = 3 periods of NREM. To address this question, we hrst introduce a 


new measure in Section 4T, the conditional power spectrum, which quantifies associations 
between power spectra and outcomes. Then, in Section |4.2[ we develop a tensor product 
model for the conditional power spectrum that extends the Bayesian spline model of Cholesky 
components of a single multivariate time series to account for dependence on both frequency 
and outcome. 

As previously mentioned, there are two approaches to conducting a Bayesian analy¬ 
sis with splines: empirical Bayes and fully Bayesian. Each approach has strengths and 
weaknesses. In the empirical Bayes approach, smoothing parameters are estimated through 
a data-driven procedure. Estimates conditional on smoothing parameters can be quickly 
computed through Fisher’s scoring or Newton-Raphson and conditional inference on the 
modeled functions (Cholesky components in our setting) can be conducted through approxi¬ 


mate “Bayesian confidence intervals” (Gu, 1992). In the fully Bayesian approach, smoothing 
parameters are treated as random variables with noninformative priors and MCMC tech¬ 
niques are used to sample from the joint distribution of coefficients and smoothing parameters 
( Speckman and Sun| 2003; Crainiceanu et ah, 2005). The sample simulated from the pos¬ 
terior distribution using MCMC provides a natural means of conducting inference on any 
function of the spectrum averaged over the distribution of the smoothing parameters, which 
accounts for uncertainty in the smoothing parameters when conducting inference. As will 
be illustrated in Section inference on squared coherence, univariate spectra, and integral 
functions thereof (all of which are nonlinear functions of modeled Cholesky components) are 
of direct scientihc interest. We develop the proposed methodology under a fully Bayesian 


framework, presenting prior distributions in Section 4.3 and the sampling scheme in Section 

1121 
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4.1 Conditional Power Spectrum 

Without loss of generality, we formulate the methodology assuming that Ui is scaled to take 
values within [0,1]. To quantify the association between the power spectrum of the time 
series Xjt and the outcome Uj, we dehne the conditional power spectrum 

OO 

/(a;, u)= Y, Gov (X,*, \ Uj = u) 6““, a; G M, n G [0,1], 

/^ = —OO 

As with the power spectrum of a single multivariate time series, the spectral matrices f{oj, u) 
are positive-dehnite P x P Hermitian matrices, and f{-,u) is a periodic and Hermitian 
function of frequency for hxed u. In a traditional spectral analysis without a cross-sectional 
variable, spectral measures such as fpq and = \fpqf / {fppfqq) are curves as functions of 
frequency. In the conditional setting, these are surfaces as functions of both frequency and 
the variable u. How these functions change with respect to u provides information as to how 
spectral measures are associated with the variable. 


4.2 Bayesian Tensor-Product Model of Cholesky Components 


As in the classical setting discussed in Section 3.2.2, where a multivariate time series observed 
from a single subject, to preserve positive dehniteness, we model the Cholesky components. 
Let 


be the modihed Cholesky decomposition of the conditional power spectrum. We use Bayesian 
tensor product models for the P^-unique Cholesky components, which decompose the bi¬ 
variate functions into products of univariate functions of u and of u. 


Bayesian models for Cholesky components as functions of u were discussed in Section 


3.2.2 Similarly, a low-rank approximate Bayesian smoothing spline model for a function of 


outcomes at the observed values can be formulated. Since the domain of the outcome values 
is [0,1], as opposed the domain of the frequency values [0,1/2], we consider the kernel 

H{ui, Uj) = {ui- v)^ {uj - v) dv. 
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Letting u = {ui,... iUn)'■, a low-rank model for functions of the outcome evaluated at the 
observed values is: 

Lna -I- (5) 

where Lh = | wj, Qh is the N x uh matrix of the hrst uh columns of H = 

VhD]^‘^VIj is the spectral decomposition of the N x N matrix H = {H{ui,Uj)}, and a ~ 
A^(0, (Jah) is independent of 6 ~ -^(0, T‘^In„)- 

To write the tensor-product model at the observed frequency-outcome points, concatenate 
components across frequency and outcome to dehne the A^M-vectors 

0k£ = Ui), . . . , Qke{l^M,Ui)} , ■ ■ ■ , {©^(^1, Mat), . . . , "Wv)}] 

forA:>£ = l,...,P — 1. Similarly dehne log^^^^ for A; = 1,..., P. The real and imaginary 
parts of Oke, and log^^^^ can then be expressed as tensor products of the spline models for 
functions of frequency (given in Equations (|^ - (|^) and outcome (given in Equation ([^) 

ke] = {Lh ^ Lj} ttrke + {Qh Lj}brke + {Lh ^ Qj} Crki + {Qh Qj}drk£ 

= {Lh ^ Lj} aik£ + {Qh Lj}bik£ + {Lh <S) Qj} Cik£ + {Qh ^ Qj}dik£ 

logV'fcfc = {Lh Lj} adkk + {Qh Lj}bdkk + {Lh Qj} Cdkk + {Qh Qj}ddkk- 

This model decomposes conditional Cholesky components into combinations of products of 
univariate functions of frequency and univariate functions of outcome. The parameters a 
are coefficients for functions that are products of linear functions of both u and u, b are 
coefficients for functions that are products of linear functions of u and nonlinear functions 
of u, c are coefficients for functions that are products of nonlinear functions of uj and linear 
functions of m, and d are coefficients for functions that are products of nonlinear functions 
of u and of u. 

4.3 Prior Distributions 

We dehne two types of prior distributions: prior distributions on coefficients conditional on 
smoothing parameters and prior distributions on smoothing parameters. The tensor product 
model naturally enables the formulation of prior distributions that regularize its components 
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as univariate functions of frequency and outcome. Letting 77 ^^^ = 6 ^^^) ^rkiy 1 Viki = 

c-M) <H)'and ri^^k = io,'dkk^Kkk,(^dkk,ddkky^ conditional on smoothing parameters, 
we assume the diagonal Gaussian smoothing priors 


Tfrke ~ ^(0> Drki) where Drki = diag((T^l 4 , 
"Hiki ~ ^(0) Diki) where = diag(cr^l 4 , 
Vdkk ~ Ddkk) where D^kk = diag(a^l 4 , 
where 1„ is the n-vector of ones. 


I3rkl^2n u 1 '^^rkl^2nji '^Srkt^njjXrij) i 

2 - 1 / 2 - 1 / 2 - 1 / 'i 

'l3iki^2nH'’ '^ikl^2nj'> 'Siki^nn^nj) 

2 - 1 / 2 - 1 / 2 - 1 / \ 

'l3dkk^2nH ^ ''ydkk^2nj'! 'Sdkk^nHX'nj) 


Prior distributions on the smoothing parameters are placed by assuming that T-yrki, 
'^5rk£i '^f)ikii '^■yikij Tsikit k > £ 1 ,..., P* 1 , T^dkki 'T'sdkki k 1 ,..., are inde¬ 

pendent Half-f(i 2 , G) random variables with pdf p{x) oc [1 - 1 - {x/GY a; > 0, 


where the hyperparameters v and G are assumed known (Gelman, 2006). Gomputationally, 


it is convenient to utilize the following scale mixture representation (Wand et ah, 2012): 

I g) ~ ^G{u/2,u/g), g ~ /G(l/2,1/G^), where IG{a,b), is the inverse Gamma dis¬ 
tribution with pdf p{x) oc exp(— 6 /a:), x > 0. The larger the value of G, the less 

informative the prior, and we set G to a large hxed value. We found analyses to be insen¬ 
sitive to the choice of G, with G = 10 and G = 10^ giving indistinguishable results in both 
simulations and in the analysis of the AgeWise data. The hyperparmeter cr^, which is the 
prior variance of the coefficients of the linear terms, is assumed to be a known large value. 
In our computations, cx^ = 10^ and cx^ = 10^ gave indistinguishable results. 


4.4 Whittle Likelihood, Sample Scheme and Inference 

Given observed time series, we dehne the DFT for the jth subject at frequency Um as 

n 

y jm = exp(-27ria;^t). 

t=i 

For large n, conditional on Uj, Yjm are approximately independent mean-zero complex Gaus¬ 
sian random variables. This provides the conditional Whittle likelihood 

N M 

L{y I /) ~ n W\f~^^^rn,Uj)\exp{-Y*^f-^{uJm,Uj)Yjm). 

j=l m=l 
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There are [{uh + 2) (nj + 2) + 3] parameters in the model of /: {uh + 2) {rij + 2) re¬ 
gression coefficients and 3 smoothing parameters for each of the Cholesky components. 
We develop a sampling scheme to sample from the joint posterior distribution of the co¬ 
efficients 77 ’s and smoothing parameters r^’s conditional on the DFT Y and the observed 
outcomes u. To aid in developing this sampling scheme, it is advantageous to consider a 
more compact notation by dehning 

Q = (^Lh ® Lj I Qh ® Lj I Lh ® Qj I Qh ® Qj^ 

so that 

= QVrki^ ^{Oke} = QViki and = QVdkk- 

Each iteration of the sampling scheme consists of three steps. First, the coefficients corre¬ 
sponding to the real and imaginary components of 0 and are sequentially sampled 
as Gaussian random variables from their conditional posterior distributions conditional on 
the current values of all other parameters. In the second step, the coefficients correspond¬ 
ing to the diagonal elements of iViikk) are drawn. The log of the conditional posterior 
distribution of these coefficients is given by 

N M ^ 

^ogp{vdkk I Vk,Ddkk) = imJtdkk ^ i^dkk^ dkk'Hdkki ( 6 ) 

j=l m=l 

where is the row of Q corresponding to the jth subject and mth frequency, Vk is a vector 
with components Vkjm depending on Y and on other parameters held hxed (its exact form 
is given in Appendix B), and = denotes equality up to a constant. Since this is not a known 
distribution, rj^kk are drawn in a Metropolis-Bastings step. The last step samples smoothing 
parameters from their posterior distributions conditional on other parameters. For ease of 
notation, in what follows we describe the sampling scheme for the case P = 3. Further details 
are given in Appendix B. After initializing all the parameters, the sth iteration, 1 < s < S, 
of the Gibbs sampler consists of the following steps. 

1. Sample the coefficients corresponding to 0: 
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(a) Draw 


~ iV(C,21.Sr2l) 

(>f.'2l|r,^u'''".«3r‘’.A'2.'‘’) ~ ^V(^a,Ei2l) 


and update O 21 

= QVr2i + iQvSi- 


(b) Draw 

Vlr31 1 r 5 V^ll > *'21 ) ^r31 

) ~ N {n^2,li^r3l) 




and update 

= Q^',% + iQVav 


(c) Draw 

U/r32 1 r ) ^22 5 -^2^32 ) 

~ ^(/^r32) Dr32) 


Vli32 \ r ) V^22 ) ^i32 ) 

~ ^ { 1 ^ 1321 ^ 132 ) 

and update 6>^J = Qrti% + ^QVill- 



The exact forms of the conditional means and covariances, ficki and T,cke, c = r,i, are 
given in Appendix B. 


2. Sample coefficients corresponding to T 
for /c = 1, 2,3 do 


Draw 


Y^dkk), where fj^kk is the maximizer of 0 and Y^dkk is 
the inverse of the observed information matrix at ffdkk- 


(b) Compute 

_ Ptffcfc I Vk,Ddkk)fT{r}'!kk^) 

P[v^lkk^ \Vk:Ddkk)fT{Tt!lk)' 

where fx is the density of the tyijfdkki ^dkk) distribution. 

(c) With probability min(l,r(")) accept otherwise rj^^^k = Vdkk^^- 

(d) Update = exp(g77^Jj. 

end 


3. Sample smoothing parameters: 
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for £ = 1, 2 do 
for /c = 2, 3 do 
Draw 

^Iril ~ ^G^(K + ^)/2, + 

~ ^G^(K + ^^)/2, c2cS/2 + z//^^;"/^) 

T-^rfc? ~ ^<^((^d + '^)/2, <^i142/2 + ^^/c/iy^) 

~ ^G{{i' + l)/2, + 1/G^) 

S'-^riz ~ + l)/2, + 1/^^) 

fi'lrL ~ + l)/2, ^/'hXl + 

end 

end 

The smoothing parameters for the imaginary and diagonal components are similarly 
drawn from inverse gamma distributions. 

Point Estimates and Credible Intervals 

The sample generated via MCMC methods provides a means of obtaining point estimates 
and credible intervals for any function of the spectrum averaged over the distribution of 
smoothing parameters through the sample mean and percentiles of the empirical distribution 
of the function evaluated at each iteration of the sampling algorithm. For instance, a measure 
of interest in the analysis of HRV is the log-spectrum from the pth period of NREM, log fpp. 
Consider 

as the estimated spectral matrix at the sth iteration corresponding to uj and Um with pth 
diagonal element f^\um,Uj). The matrix Q^'^\um,Uj) has kith element 

^kljm = QjmViM + i QjmViW k > i = 1, . . . , P - 1, 

and Uj) has kkth element exp (^g'kkVdkk^ ■ If ^ iterations of the sampling algorithm 

are run with a burn-in of Sq, then an estimate of log fpp{um, Ui) can be computed as the 
mean of the values {log/p? {um,Ui ); 5*0 < s < S'! , and a 95% credible interval computed 
as their 2.5 and 97.5 empirical percentiles. 
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Scientific interest also lies in measures collapsed across frequency. For example, as will 
be discussed in the following section, in the analysis of HRV, collapsed power within the 
high-frequency band (HF) between 0.15-0.40 Hz 

/ •.40 

fpp{u,u)du (7) 

15 

provides a measure of autonomic nervous system activity during the pth NREM period 

among people with a TIB of u. Letting fp^^^\u£) = W~^ ^ i 5 <a;,„< 4 o where W 

is the number of Fourier frequencies within the HF band, an estimate of fp^{ui) is given by 

r HF fs) '1 

the mean of the values < fp ^ yUi); So < s < S >, and a 95% credible interval is given by 
their 2.5 and 97.5 empirical percentiles. 


5 Application to the AgeWise Study 


We used the proposed methodology to analyze the association between TIB and the power 
spectrum of the hrst three periods of NREM from N = 108 AgeWise subjects, as described 
in Section]^ The method was fit using uh = nj = 10 basis functions, with hyperparameters 
G = = 10^, and for S = 3000 iterations of the MCMC algorithm with a burn-in of 

So = 500. Note that, in this example, there are a total of {uh + 2) (rij + 2) = 1296 
coefficients and 3P^ = 27 smoothing parameters. The average run time per iteration was 
5.46 seconds with a standard deviation of 0.24 seconds using the program that is available 
on the journal’s website in Matlab 2016b and macOS Sierra vlO.12.1 on a 2.9 GHz Intel 
Core i7 processor with 16 GB RAM. 


Although all desired analyses are obtained from one MCMC chain, to aid the biological 
and clinical discussion of the results, we present the analysis in three stages. First, in Section 


5.1 we examine the estimated period-specihc spectra and squared coherences as frequency- 


outcome surfaces. In the subsequent two stages, we explore power collapsed within certain 
frequency bands as functions of TIB: hrst for power within periods in Section 5.2[ then 


for coherence between periods in Section |5.3[ The results from these analyses provide new 
insights into biological underpinnings of spending too little or too much time in bed. In 
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particular, our analysis suggests that (i) short TIB is connected to elevated stress and arousal 
within-periods of NREM towards the end of the night and (ii) long TIB is associated with 
a persistence in arousal in the beginning of the night. 

5.1 Analysis of the Conditional Spectrum 

Point estimates of the within-period conditional log-spectral surfaces, \og{fpp{oj,u)}, and of 
the cross-period conditional logit squared coherence surfaces, 

logit = log [plg{uj,u)/ { 1 -plg{uj,u)}] , 

are displayed in Figure These estimates are plotted on the logarithmic and logistic scales, 
respectively, to aid visualization. The conditional spectra at each of the hrst three periods of 
NREM and the squared coherence between NREM 1 and 2 display different characteristics 
within low frequencies that are less than 0.15 Hz compared to higher frequencies between 
0.15-0.40 Hz. 


From a biological perspective, these results are not surprising and produce interpretable 
measures. As was discussed in Section the autonomic nervous system is classically divided 
into two branches: the parasympathetic branch that is responsible for activities related to 
resting and digestion and the sympathetic branch that is responsible for the flight-or-£ght 
response. Researchers have shown that power within the high frequency band (HE) within 
0.15-0.40 Hz provides a measure of parasympathetic nervous system activity and that power 
within the low frequency band (LF) between 0.04-0.15 Hz is a measure of the combined 
modulation of both the sympathetic and parasympathetic nervous systems. Consequently, 
the ratio of power from low frequencies versus high frequencies (LF/HE) can be interpreted as 
a measure of sympathetic modulation relative to parasympathetic modulation. Blunted HE 
and elevated LF/HF power are often interpreted as indirect measures of physiological arousal 
and psychological stress ( Hall et 2004, 2007). To obtain inference on associations between 
these measures and TIB, in the next two subsections we examine power and coherence 
collapsed within these bands as functions of TIB. 
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NREM 1 


NREM2 


NREM 3 



NREM 1-2 


NREM 2-3 


NREM 1-3 





Figure 2: Estimated conditional log-spectra for each period of NREM (top panel) and esti¬ 
mated logit of conditional coherence between each period (bottom panel). 


5.2 Analysis of Within-Period Power 


We consider two collapsed measures of within-period power. In addition to HE previously 
defined in Equation ([^, we also consider LF/HF as 

= |y fpp{u,u)du^ ^ fpp{u,u)du^. 

Estimates and 95% pointwise credible intervals for these two measures as functions of TIB 
are displayed in Figure]^ for each period. 

HF power is relatively constant across TIB during NREM 1, while participants with a 
TIB of less than 400 minutes have decreased HF power during NREM 2 and 3 compared 
to those who spend more time in bed. Further, those who have an exceedingly small TIB 
display increased LF/HF power during NREM sleep compared to those who spend more TIB, 
especially during NREM 3. These characteristics are indicative of heightened physiological 
arousal and psychological stress. 
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NREM 1 
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Figure 3: Estimated conditional HF (top panel), fp^, and LF/HF (bottom panel), , 

as functions of TIB with 95% pointwise credible intervals for each period of NREM. 


Sleeping less than 7 hours per night has been shown to be associated with a multitude 


of negative health effects, including increased mortality (Buysse, 2014). The results of our 
analysis provide a potential pathway through which short sleep, which is inherently bounded 
by TIB, is connected to well-being: through increased stress and arousal towards the end of 
the night. 


5.3 Analysis of Cross-Period Coherence 


To investigate connections between cross-period coherence and TIB, we consider conditional 
HF band-squared coherence 


Pirin) = 


'.40 


fpq{u;,u)du 


/{/rw/rw} 


^.15 / 

and display estimates on the logit scale, logit = log (l —p^^^'^)], in the 

top panel of Figure To better understand how changes in TIB are associated with HF 
coherence, we also examine first derivatives. 


Dpg^iu) = d [plf^iu)] /du, 
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Figure 4: Estimated logit of integrated HF coherence (top panel), logit between 

each NREM period as functions of TIB and their Erst derivatives (bottom panel), 
with pointwise 95% credible intervals. 


whose estimates are displayed in the bottom panel of Figure]^ We hnd that the derivative 
of HF coherence between NREM 1 and 2 is positive for TIB greater than 500 minutes. This 
indicates that excessive increases in the amount of time spent in bed are associated with 
increased coherence in parasympathetic activity in the beginning of the night. 


The relationship between excessive TIB and ill-health led Youngstedt and Kripke (2004) 
to propose modest sleep restrictions to increase quality of life and survival, especially for 
older adults, who tend to spend more time in bed as compared to younger adults. However, 
these restrictions must be used with great care as they can potentially lead to negative 


health effects (Reynolds HI et ah, 2010 Reynold et ah, 2014). Our results demonstrate that 
excessive TIB is associated with a coherence in parasympathetic activity in the beginning 
of the night that is not present in moderate and short TIB. A possible explanation for this 
relationship is that extensive TIB can cause an increase in the amount of time spent awake 
while in bed, or lead to fragmented sleep. The roles of and relationships between physiological 
activity during different sleep cycles could change as sleep becomes more fragmented. These 
Endings provide some of the Erst potential insights into the biological pathway through which 
excessive TIB can be connected to negative health, which can potentially be used to inform 


26 

























optimal sleep restriction strategies in older adults. 


6 Final Remarks 


This article introduces a novel approach to analyzing associations between multiple time 
series and cross-sectional outcomes when data are observed from multiple subjects. A new 
measure of association, the conditional power spectrum, is introduced and its Cholesky 
components are modeled as bivariate functions of frequency and cross-sectional outcome. A 
MCMC algorithm is developed for model fitting allowing for inference on any function of the 
power spectrum. The method was motivated by a sleep study and uncovered connections 
between excessive time in bed and heightened arousal and stress that could not have been 
uncovered through traditional methods. 


We conclude this section by discussing three extensions to the proposed methodology. 
First, the model is formulated to investigate the association between power spectra and a 
single cross-sectional variable. The model could easily be extended through higher-order 
tensor product models to include multiple variables, such as the amount of time it takes to 
fall asleep and the number of awakenings during the night. However, such a model would 
provide inference on the effect of these variables on the power spectrum conditional on 
the other variables, complicating interpretation when these variables are highly correlated. 
Future work will explore an interpretable canonical correlation type dimension reduction of 
a collection of correlated variables and multivariate spectral matrices, which can be viewed 


as a multivaraite extension of Krafty and Hall (2013). Second, our application focused on 
HRV, due to the insights that it provides into autonomic nervous system activity. One could 
also explore the spectral analysis of other PSG channels, as well as the simultaneous coupling 
of channels. However, each channel of the PSG is sampled at a different rate. The second 
extension will develop conditional spectral analysis of time series with different sampling 
rates. Finally, since we were motivated by the analysis of HRV during epochs within NREM 
that are approximately stationary, we focused on stationary time series. For more highly 
sampled signals such as EEG, this assumption is not valid. A conditional time-frequency 
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analysis for signals that are locally stationary will also be explored. 


Code 

A Matlab program for implementing the proposed methodology and a hie demonstrat¬ 
ing its use are available at https://www.mathworks.com/matlabcentral/hleexchange/61186- 

mcbspec. 
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Figure A.l: Simulated conditional MA(2) epochs of length n = 300 with Uj = 0.04 and 

Uj = 1 . 

Appendices 


A Simulation Study 


To illustrate the proposed model and to investigate its empirical properties, we consider the 
P = 3 dimensional second order moving average model ( MA(2) ) 


~ A J — 1; ■ ■ ■ 5 A^5 t — 1,... ,n, 


where 0i = —I, 02 = 0.6/, / is the 3x3 identity matrix, and ejt are independent N [0, (uj)] 

random variables with 


n (u) = cr‘^{u) 


1 p{u) p{u) 
p{u) 1 p{u) 
pin) p{u) 1 


= (2 — and p{u) = 0.6 + 0.25 cos (ttu). Two simulated epochs of length n = 300, 
one with Uj = 0.04 and one with Uj = 1, are displayed in Figure A.l The conditional 
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Figure A.2: Conditional log-spectra from the MA(2) model (top panels) and estimated 
conditional log-spectra from a random sample of time series of length n = 300 from N = 25 
subjects (lower panels). 


spectrum is given by 

f{oj, u) = 0(a;)ff (u) 

where ©(cn) = I + ©i exp(—27ria;) -f ©2 exp(—dvrio;). The log-spectra, \og[fpp{oj,u)], and 
their estimates under the proposed procedure from a random sample of A^ = 25 independent 
epochs of length n = 300 are displayed in Figure |Al2| Plots of the logit-squared coherence, 
logit [ppg{uj, u)~\, and their estimates are displayed in Figure A.3 The band-collapsed mea¬ 
sures and along with their estimates and 95% credible intervals. 


are displayed in Figure A.4 


We simulated 100 random samples of conditional MA(2) time series of length n from N 
subjects with Uj = j/N for the four combinations of n = 300,500 and N = 25,50. The 
estimation procedure was run using nj = 10, uh = 5, and for 2000 iterations of the MCMC 
algorithm with burn-in of 500 iterations. To investigate the sensitivity of the proposed 
estimation procedure with respect to hyperparameters, the procedure was run twice: for 
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Figure A.3: Conditional squared coherence from the MA(2) model (top panels) and estimated 
conditional squared coherence from a random sample of time series of length n = 300 from 
N = 25 subjects (lower panels). 



n = 300 

n = 500 

AM 


0^5 

0^2 



(0.02) 

(0.03) 

N 

= 50 

1.29 

3.40 



(0.03) 

(0.28) 


Table A.l: Mean (standard deviation) run time per iteration of the sampling algorithm in 
seconds. 


G = 10® and for G = 10^°. Table A.l reports the mean and standard deviation of run times 
per iteration for each setting using the program that is available on the journal’s website in 


Matlab 2016b and macOS Sierra vlO.12.1 on a 2.9 GHz Intel Core i7 processor with 16 GB 


RAM. 


To investigate the performance of the proposed procedure for conducting inference on 
band-collapsed measures as functions of outcome, we computed pointwise 95% credible in¬ 
tervals for the nine band-collapsed curves ■, ^ 

P 23 ^^, and The mean and standard deviation of pointwise coverage probabilities in¬ 


tegrated across u are given in Table A.2 The integrated coverage was near the nominal 
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Figure A.4: HF band power (top panels), LF/HF band power (middle panels), and HF band 

squared coherence (bottom panel) (—), along with point estimates ( - ) and 95% credible 

intervals (- ■ -) from a random sample of = 25 conditional MA(2) time series of length 
n = 300. 
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n 

N 

G 

fr 

fHF 

J2 

fHF 

i.S 

rLF/HF 

rLF/HF 

12 

rLF/HF 

h 

2,HF 

Pl2 

2,HF 

P2^ 

2,HF 

Pi?, 

300 

25 

10^ 

.943 

(.112) 

.968 

(.073) 

.958 

(.097) 

.969 

(.084) 

.965 

(.078) 

.961 

(.093) 

.968 

(.058) 

.951 

(.083) 

.963 

(.078) 



10^0 

.943 

(.112) 

.968 

(.073) 

.958 

(.097) 

.969 

(.084) 

.965 

(.078) 

.961 

(.093) 

.968 

(.058) 

.951 

(.083) 

.963 

(.078) 

500 

25 

10^ 

.968 

(.070) 

.949 

(.087) 

.969 

(.070) 

.965 

(.097) 

.974 

(.062) 

.967 

(.069) 

.968 

(.068) 

.950 

(.086) 

.963 

(.083) 



10^0 

.968 

(.070) 

.949 

(.087) 

.969 

(.070) 

.965 

(.097) 

.974 

(.062) 

.967 

(.069) 

.968 

(.068) 

.950 

(.086) 

.963 

(.083) 

300 

50 

10^ 

.948 

(.107) 

.961 

(.077) 

.960 

(.073) 

.962 

(.092) 

.945 

(.110) 

.945 

(.106) 

.962 

(.082) 

.961 

(.073) 

.960 

(.081) 



10^0 

.948 

(.107) 

.961 

(.077) 

.960 

(.073) 

.962 

(.092) 

.945 

(.110) 

.945 

(.106) 

.962 

(.082) 

.961 

(.073) 

.960 

(.081) 

500 

50 

10^ 

.967 

(.070) 

.956 

(.079) 

.971 

(.06^) 

.954 

(.102) 

.955 

(.093) 

.963 

(.078) 

.959 

(.080) 

.953 

(.088) 

.969 

(.059) 



10^ 

.967 

(.070) 

.956 

(.079) 

.971 

(.06J,) 

.954 

(.102) 

.955 

(.093) 

.963 

(.078) 

.959 

(.080) 

.953 

(.088) 

.969 

(.059) 


Table A.2: Mean (standard deviation) coverage of 95% credible intervals for band-collapsed 
measnres from 100 random samples of N conditional MA(2) time series of length n nsing 
hyperparamter G. 


95% level for each component, ranging between 94.3%-97.4%. Coverage probabilities nnder 
different tnning parameters were indistingnishable. 

To compare the performance of the proposed procednre to existing approaches, we also 
computed two two-stage estimators of within-period band-collapsed measures. In the hrst 
stage, periodograms were calculated and summed within HF and LF bands for each of the 
P = 3 series for each of the N subjects to obtain raw subject-specihc estimates. The raw 
subject-specihc HF and LF/HF estimates were then smoothed across u. For the hrst estima¬ 
tor, smoothing was achieved by htting a cubic smoothing spline with smoothing parameter 
selected through generalized cross-validation (GCV) (|Gu 2013). For the second estimator. 


smoothing was achieved through local linear regression with plug-in bandwidth (Loader 


1999). We dehned the integrated square error (ISE) of an estimate of as 


- /rw 


du. 
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The ISEs for and Ppq^^ were similarly defined. The mean and standard deviation 


of the ISEs are presented in Table |A.3[ As expected, the ISE of each estimator improved 
with an increase in either n or N. The insensitivity of the proposed procedure to choice of 
hyperparameter that was observed through indistinguishable coverage probabilities was also 
observed in the ISE; the ISE under G = 10® and G = 10^*^ were identical up to at least three 
signihcant digits. In each setting, the proposed estimator had smaller mean ISE compared 
to the two-stage procedures. 


37 



n 

N 

Estimator 

fr 

fHF 

12 

fHF 

is 

nLF/HF 

Jl 

rLF/HF 

12 

pLF/HF 

Js 

300 

25 

Bayes: 10^ 

6.76 

(4.96) 

6.75 

(6.03) 

6.57 

(6.30) 

2.98 

(2.23) 

3.28 

(2.45) 

4.31 

(3.34) 



Bayes: 10^° 

6.76 

(4.96) 

6.75 

(6.03) 

6.57 

(6.30) 

2.98 

(2.23) 

3.28 

(2.45) 

4.31 

(3.34) 



2-Stage: Spline 

12.11 

(12.70) 

13.14 

(18.56) 

12.82 

(12.20) 

8.75 

(8.29) 

10.06 

(8.29) 

8.92 

(7.00) 



2-Stage: LOESS 

10.31 

(6.97) 

11.25 

(14.06) 

11.68 

(7.89) 

10.01 

(7.36) 

10.86 

(7.36) 

10.18 

(7.97) 

500 

25 

Bayes: 10^ 

3.75 

(2.94) 

4.57 

(3.54) 

3.75 

(2.82) 

1.86 

(1.68) 

1.79 

(1.29) 

2.62 

(1.89) 



Bayes: 10^° 

3.75 

(2.94) 

4.57 

(3.54) 

3.75 

(2.82) 

1.86 

(1.68) 

1.79 

(1.29) 

2.62 

(1.89) 



2-Stage: Spline 

7.87 

(7.76) 

8.05 

(8.63) 

7.15 

(7.27) 

5.26 

(4.35) 

5.04 

(4.80) 

4.60 

(3.44) 



2-Stage: LOESS 

7.93 

(6.32) 

7.69 

(5.29) 

7.21 

(5.61) 

5.94 

(3.84) 

5.34 

(3.50) 

5.42 

(2.87) 

300 

50 

Bayes: 10^ 

3.33 

(2.69) 

3.61 

(2.84) 

3.59 

(2.97) 

1.51 

(1.37) 

1.83 

(1.39) 

2.51 

(1.81) 



Bayes: 10^'^ 

3.33 

(2.69) 

3.61 

(2.84) 

3.59 

(2.97) 

1.51 

(1.37) 

1.83 

(1.39) 

2.51 

(1.81) 



2-Stage: Spline 

6.64 

(7.70) 

7.09 

(10.34) 

6.84 

(6.63) 

4.88 

(4.15) 

5.39 

(4-41) 

5.57 

(4.35) 



2-Stage: LOESS 

5.50 

(4.00) 

5.60 

(6.35) 

5.91 

(4.21) 

5.31 

(3.64) 

5.62 

(3.48) 

5.79 

(3.85) 

500 

50 

Bayes: 10^ 

1.74 

(1.25) 

2.12 

(1.50) 

2.06 

(1.36) 

1.00 

(0.83) 

1.13 

(1.13) 

1.44 

(0.92) 



Bayes: 10^'^ 

1.74 

(1.25) 

2.12 

(1.50) 

2.06 

(1.36) 

1.00 

(0.83) 

1.13 

(1.13) 

1.44 

(0.92) 



2-Stage: Spline 

4.52 

(5.00) 

4.87 

(5.31) 

3.84 

(4.31) 

2.87 

(2.56) 

3.19 

(3.07) 

2.84 

(2.32) 



2-Stage: LOESS 

3.85 

(3.13) 

3.89 

(2.89) 

3.47 

(2.41) 

2.81 

(1.92) 

3.04 

(2.18) 

2.79 

(1.57) 


Table A.3: Mean (standard deviation) of the integrated square error (ISE) of band-collapsed 
measures from 100 random samples of N independent conditional MA(2) time series of length 
n. Estimates were obtained using the proposed procedure with tuning parameter G = 10^ 
(Bayes: 10^) and G = 10^*^ (Bayes: 10^°) and two-stage estimators using smoothing splines 
(2-Stage: Spline) and local linear regression (2-Stage: LOESS). Values are reported xlO^ 
for HE measures and xlO^ for LF/HF measures. 
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B Details of the Sampling Scheme 


In this appendix we provide more details about the sampling scheme outlined in Section 4.4, 
assuming P = 3. As in Section 4.4, is the row of Q corresponding to uj and Um- The 
DFT of the pth series from the jth subject at Fourier frequency ojm is denoted by Ypjm- The 
kl element of Q{um,Uj) defined in Section 4.2 is expressed as 

Okljm = q'jmVrke + * Q'jmVike, k > i = 1,..., P - 1, (B.l) 


where the i in the second term on the right-hand side of (B.l) is the unit imaginary number. 
The diagonal elements of are expressed as ~ ^^viq'jmVdkk)^ ^ = 1,... ,P. 

To aid presentation and simplify notation, the superscript for iteration number is suppressed 
and all derived distributions are conditional on the current values of all other parameters. 


Drawing the Basis Function Coefficient Vectors 


The conditional posterior distribution of c = r, i, fc > £ = 1,..., P — 1, is multivariate 
normal, T^cki)- In what follows we provide expressions for and Scm- 

Mean vectors and covariance matrices for 77^21 Vai 

N M 

^c21 ~ 2 + -^^21! C = r,i 

j=l m=l 
N M 




2lt^r21 




n* V 



j=l m=l 
N M 




jm- 


j=l m=l 
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Note that ipiijm ^sijm are evaluated at their current values. 
Mean vectors and covariance matrices for and 

N M 


y-1 

^c31 


^rSl/^rSl 


+-^csu c — r,i 

j=l m=l 
N M 




-e 


* 

21 jm 


Y, 


2jm 



j=l m=l 
N M 




11 jm 

j=l m=l 






jm- 


Mean vectors and covariance matrices for 77^32 and 77^32 

N M 

^c32 = ^ ' 1 ^ 227 'ml^ 3 jm| + -^£32) C = r,i 

j=l m=l 

N M 

^r32/^r-32 = 2EE '^22jm 

jf=l m=l 

N M 

^i32t^i32 — 2EE '^22^jm 

jf=l m=l 


The basis function coefficient vectors k = 1,..., P, are drawn from p( 77 ^^;;. | Q, Vk, Ddkk ), 
given in Section 4.3. The entries v^jm of Vfc for /c = 1, 2, 3 are as follows. 


I^ljm 


^2jm 

'^3jm 


\Yljm\‘^ + \O 2 ij m Y2jm\ T 


2 


|h2jmP + \d32jmY3jm\^ — 23?1032jmh2jmf^3im} • 


1^3,mP. 


The vectors rj^kk: ^ = 1,... ,P, are generated independently via a Metropolis-Hastings step 
with a multivariate t proposal distribution, tyifidkki ^dkk), where 

Vdkk = arg max logp{‘ndkk I Q,Vk,Ddkk) 

Vdkk 


and 


^dkk — 




dVdkkdVa 


log p ( 77 , 


dkk 


dkk 


Q1 ^ki 


-1 

dkk dkk 
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The gradient and Hessian of \ogp{r}^i.f^ \ Q,Vk, Ddkk) are given by 


and 


respectively. 


N M 


Y.Y. 1 - Vkjm exp{q'j^rtakk)\ 9jm ^ dkk^dkk 


j=l m=l 


N M 

-EE 

j=l m=l 


Drawing the Smoothing Parameters 

Details are given below for the smoothing parameters associated with the real part of 
6ke{(^,u). The details for the rest of the smoothing parameters are similar. The smoothing 
parameters r'^^ke^ '^^rki '^Srki drawn independently for A; > £ = 1,..., P — 1, as follows. 

^Irfez ~ IG{{nb + z/)/2, Kkibrki/2 + v/g^rki) 

T^rkl ~ IG{{nc + u)/2, C^^iCrkl/2 + u/g^rkl) 

Tsrki IG{{nd + iy)l2,d’^,^idrki/2 + iy/gsrki), 

where 

ggrki - IG{{i. + 1)/2,p/tI,, + 1/G^) 
g,^ki - /G((z. + 1 )/2 ,z./<,, + 1 /G2) 

gsrki ~ IG{{u + l)/2,u/rl,, + l/G^). 
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